technical report section 

Laval POSTGRADUATE SCHOOL 

MONTEREY. CALIFORNIA 93940 



NPS-58G172071A 



NAVAL POSTGRADUATE SCHOOL 

Monterey, California 




THE DEVELOPMENT OF A HOMOGENEOUS NUMERICAL 
OCEAN MODEL FOR THE ARCTIC OCEAN 

by 

J. A. Galt 
July 1972 

Approved for public release; distribution unlimited 



FEDDOCS 

D 208.14/2:NPS-58GL72071A 




NAVAL POSTGRADUATE SCHOOL 
Monterey, California 



Rear Admiral A. S. Goodfellow, USN Milton U . Clauser 

Superintendent 



ABSTRACT: 

A numerical ocean model driven by surface stress and a source- sink 
distribution is developed for a homogeneous ocean. Non-linearities, 
lateral friction and bottom friction are included. The basin shape can be 
varied to accommodate a large variety of configurations. Variable bathy- 
metry and sources/sinks around the perimeter are included. The numerical 
scheme is conditionally stable and has second order accuracy in space 
and time . 

A number of test cases are run to explore the dynamic significances 
of the various processes represented. The possible influence of these 
processes on the circulation of the Arctic ocean are discussed. 



TABLE OF CONTENTS 



Section I - Introduction and Report Summary 4 

Section II - Development of the First Stage Numerical Model . 8 

Section III - Preliminary Results 21 

Section IV - Future Model Plans 32 

References 34 

Appendix - Use of Computer Program 36 



LIST OF FIGURES 



Figure 1 



Figure 2 
Figure 3 
Figure 4 
Figure 5 
Figure 6 



Figure 7 



. Chart cf the Arctic Basin with the 100, 1000 and 2000 fathom 
depth contours drawn. Cross indicates location of the North 
Pole . 

. Grid pattern used in finite difference scheme and numbering 
system used for calculational molecule. 

. Configuration of the model used to check out the relaxation 
of the stream function. 

. Streamlines of source-sink driven flow with zero vorticity in 
the interior and uniform depth 

. Streamlines of source-sink driven flow with uniform negative 
vorticity in the interior and uniform depth . 

. Streamlines of flow driven by uniform stress curl in an irregu- 
lar shaped ocean with a central ridge. (The ridge runs from 
top to bottom in the figure.) 

. Streamlines of flow driven by a simple source-sink distribution 
for an irregular shaped ocean with a central ridge. (The ridge 
runs from the top to the bottom in the figure . ) 



3 



Section I - Introduction and Report Summary 

The general equations that describe the flow of both the atmosphere 
and the ocean are well known and are available to geophysicists for the 
investigation of a large variety of circulation problems. For many cases 
of interest analytic solutions to these general equations are not possible 
because of significant non-linearities in the equations, or because of 
irregular geophysical boundary configurations, or both. Although analytic 
solutions are not generally available it is often possible to obtain use- 
ful approximate solutions using numerical techniques. This numerical 
modeling has become a powerful tool for the investigation of both the 
ocean and the atmosphere. 

Large scale numerical models available for the description of many 
features of the atmosphere and large portions of the oceans can be divided 
into roughly two groups. The first group is exploratory in nature and is 
used to investigate the significance of various physical processes. The 
prime emphasis is to understand just why a particular system responds as 
it does and what effects variations in forcing have on the outcome. The 
results from this type modeling may, or may not, be physically realistic 
but in all cases should elucidate the characteristics of the flow to be 
expected, its sensitivity to various inputs and significant correlations 
that are likely to exist. The second group of numerical models is generally 
more advanced and is designed for forecasting geophysical fluid flow on 
a more or less real time basis. The relative stage of refinement and the 
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large amounts of high quality data required for initial conditions has in- 
hibited the development of these type models for the ocean and essen- 
tially all of the forecasting models routinely used deal with the atmosphere. 

The object of this research was to begin a numerical exploration on 
the large scale circulation of the Arctic Ocean. 

In the past numerical models of the Arctic Ocean have concentrated 
their attention on the sea ice that overlies most of the Ocean (Campbell, 
1965) and tried to answer questions concerning the drift and climatological 
permanence of the pack ice (Maykut and Untersteiner, 1971). With the 
exception of Campbell's use of a simplified ocean model under an ice- 
layer model no significant effort has been directed towards a numerical 
study of the actual flow of the Arctic Ocean's waters. 

At the onset of this research it must be admitted that very little 
is actually known about Arctic Ocean dynamics and that difficulties can 
be anticipated in deciding on appropriate boundary conditions and input 
parameters for the model. In many cases field data from the Arctic is 
lacking, or inadequate to give the required stress fields or inflow-outflow 
conditions with sufficient accuracy. This then demands that best esti- 
mates of boundary conditions serve as a tentative guide and that wide 
ranges of parametric inputs actually be investigated. The resulting 
numerical exploration yields considerable insight into the relative signi- 
ficance of various oceanographic parameters. 

In addition to the uncertainties related to the lack of actual input 
data the development of a new numerical model has potential difficulties 
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inherent in the finite difference mathematic's that must be used. Both 
these problems can best be addressed by the careful development of the 
model from simple to more complex cases in a stepwise progression and 
with each step being checked against any available data (from the field 
or analytic considerations). This model is carried through this sort of 
genesis. 

The first step in the model to be used for the Arctic assumes homo- 
genous water and variable depth which is in some way similar to a model 
used by Holland (197 6). The grid system used is based on a triangular 
plan similar to some systems used by Williamson (1968) and Sadourny, 
Arakawa and Mintz (1967). Both lateral and bottom friction are included 
in the model. The flow is driven by stress applied at the surface (simu- 
lating the wind or ice stress) and by source-sink distributions around 
the edge (that simulates major channels between the Arctic and other 
oceans). The details of the development for this first step in the model 
are given in the next section. The initial check out and experimentation 
with this first stage in the development of an Arctic Ocean circulation 
model has led to some interesting results. In particular: 1) negative 
curl introduced into the flow results in current patterns resembling the 
Beaufort Gyre; and 2) the Lomonosov Ridge (Fig. 1) acts like a dynamic 
block which may greatly increase the significance of the circulation caused 
by the source-sink distribution. There is some indication from field data 
reported by Muench (1970) and Thorndike (1971) that the processes indi- 
cated in the model have some counterpart in actual circulation observed 
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Figure 1. Chart of the Arctic Basin with the 100, 1000 and 2000 fathom 
depth contours drawn. Cross indicates location of the North 
Pole. 



in the Arctic Ocean. A more detailed discussion of these results are 
presented in section three of this report. 

Future extension and development of the model will be continued in 
stages. In the immediate future the present model will run with greatly 
increased resolution to delineate more of the complicated bathymetry of 
the Arctic Basin. After those results have been carefully explored the 
stratification appropriate to the Arctic Ocean must be introduced into the 
modeling. The results for the stratified ocean model will indicate how 
the next step (the inclusion of more complicate thermodynamic exchanges) 
can best be approached. A more detailed discussion of future modeling 
plans is given in section four of this report. 

Section II - Development of the First Stage Numerical Model 

When one considers the actual complexity of the circulation in the 
Arctic Ocean and the wide variety of numerical modeling techniques 
available for its study it is by no means obvious which path will lead to 
maximum returns. Each of the presently used numerical models has its 
advantages and limitations. In an attempt to address this question in a 
rational way it was decided to start with the simplest numerical model that 
could simulate what are thought to be the dominant forcing and geo- 
morphology in the Arctic. 

There can be little doubt that much of the large scale circulation of 
the Arctic Ocean is wind driven either directly, or indirectly through an 
ice cover that acts as some sort of coupling element (Campbell, 1965). 
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The details of how the ice cover couples the atmosphere to the ocean and 
to what extent it filters the time and space variations are unknown. There 
is however, a substantial effort directed towards this problem (Untersteiner, 
and Fletcher, 1971) and some progress can be anticipated. For the pre- 
sent model these interesting questions can not be treated in detail and 
thus the stress on the top of the ocean will be considered a known field, 
externally specified. 

Studies by Coachman and Barnes (1961, 1963), Aagaard (1966) all 

indicate that there may be dynamically significant exchange between the 

Arctic Basins and the adjacent portions of the world ocean. For this 

\ 

reason the numerical model incoporates sources and sinks of water around 
its perimeter to simulate the major channels between the Arctic and the 
adjacent parts of the ocean. 

Looking at Figure 1 . , it is seen that a large fraction of the Arctic is 
covered by the relative shallow Siberian Shelf, while the deeper portion 
of the Arctic Ocean is divided into two basins (both over 4000 meters 
deep) by the Lomonosov Ridge. In an attempt to retain some of the dynami- 
cal effects caused by these large bathymetric variations the model includes 
variable depth. 

It is likely that the density variations in the Arctic Ocean effect the 
flow. For example a full treatment of the movement in the Atlantic layer 
(Coachman and Barnes 1963) will certainly require a baroclinic model. On 
the other hand there are some actual current measurements from the Arctic 
(Nikitin and Demyanov, 1965) (Galt, 1967) (Coachman, 1969) that indicate 
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a substantial barotropic, or depth independent component to the currents. 
This coupled with a consideration of the great increase in complexity re- 
quired for variable density models suggests that for the initial numerical 
exploration a homogeneous or barotropic formulation be used. 

The circulation in the Arctic is not well enough known to come up 
with an accurate appraisal of the significance of frictional forces. It 
seems likely that near source-sink points lateral friction could effect the 
flow. Over the large area covered by the Siberian Shelf it is quite possible 
that bottom friction might also be significant. Accordingly both lateral 
and bottom friction were included in the model and it was anticipated 
that some range of frictional parameters would be investigated to test 
for significance. 

To develop a model with the characteristics described above the 
following integrated form of the equations of motion are used: 



Du f ldP2 Rut 

— - fv = - — - — +Kv u - — + — 
Dt p 3x h ph 



( 1 ) 



Dv f 1 3P . „ 2 Rv T* 

— + fu = - - — + Kv v - — + — 
Dt p dy h ph 



( 2 ) 



Where the dependent variables u and v are the horizontal components of 

velocity which are independent of depth. The density, p, is a constant. 

r and r , the components of the wind stress, h, the depth, and f, the 
x y 

coriolis parameter, are functions of position. K and R are constants 
that specify the effectiveness of the horizontal and vertical frictional 
forces respectively. 
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In addition to these equations of motion we have the continuity 



equation: 



h (hu > + k (hv) = 



A transport stream function is 



-hu-|* 

ay 



o 

introduced such that: 



( 3 ) 



(4) 



dib 

hv = T~~ 

dx 

The pressure can be eliminated from equations (1) and (2) and after some 
manipulation the following vorticity equation is obtained: 



- (v x 0ic) • v 
ar h 

= KV 2 ^ - ^ [4 + 70 

where: 

^~dx dy 

_ d_ ^ r _d_ 

^ 1 dx ] 9y 

r = r i + r j 
x y 

—* — > — > 

and i, j, k are the unit vectors in the right handed x, y, z coordinate 
system . 

From the above the following relationship between £ and 0 is obtained 

7 (^ 90 ) = £ ( 6 ) 

Equations (5) and (6) can now be solved for the vorticity and stream 
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function provided that the proper initial and boundary conditions are 
given. In particular the following must be specified: 

a) 0 - given within the region of interest at t = 0 

b) 0 - given on the boundary of the region for all time 

c) £ - given within the region of interest at t = 0 

Note that condition b) is equivalent to specifying the source-sink distri- 
bution around the edge of the model. 

Equations (5) and (6) can be non-dimensionalized by introducing the 
following new variables: 

t=t'(^) f = f'F 

x = x' A y = y ' u 

h = h'D (7) 



0 = 0'0 o 

T = T ' ( — ) 

A 

Where the primed quantities are all non-dimensional, F is the average value 
of the Coriolis parameter, A is the finite difference grid spacing, D is 
the average depth of the model and 0 q is equal to the max range of stream 
function on the boundary of the model , or the maximum value expected for 
the wind driven portion of the flow, which ever is the most convenient. 
Using the new variables defined above equations (5) and (6) become; 



- (vx >l)'k) • v( ^ r f ~ ) 
dt h 

= 8 v 2 £' - J [£'+90' • v(^,)3 + vx(£) 



( 8 ) 
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(9) 



v(j". vib') = V 

Where: 

a. ^z DF 

. _ K 

A 2 ? 

R 

y DF 

These three non-dimensional parameters govern the character of 
the solution. For example setting a = 0 removes the non-linear advective 
terms from the model. The size of /? determines how important lateral 
friction is in the solution and y scales the importance of bottom friction. 

The finite difference grid that equations (8) and (9) are solved on is 
made up of a one dimensional array of N x M points. They are arranged 
as shown in figure 2. Grid points are numbered sequentially, left to right 
starting in the top row. This means that grid point in the neighborhood 
of the point L are given as shown in figure 2. 

Two additional one dimensional arrays are used to specify the extent 
of the interior domain. These are integer arrays labeled IA and JA of 
dimension (M-2) and are defined so that a sweep of interior model points 
is obtained via the following algorithm. 

DO 2 K=1 , M-2 

I = IA(K) 

I = JA(K) 

DO 1 L = I, J 

1 Statement on Interior Point (L) 



2 Continue 
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used for calculations! molecule . 



The depth in the model is specified in terms of an average depth and 
an array of deviations from this average. The amount of bathymetry to 
be used in any particular run is specified by a depth factor which is in- 
put at run time. The following depth variables are thus defined. 

HM = average water depth 
HFAC = depth factor 

H (DIMENSIONED TO N x M) = array of depth deviations. 

The actual depth used in any particular run is calculated with the following 
algorithm 

DO 3 I = 1 , NM 
3 H(I) = HM + HFAC * H(T) 

Obviously HFAC = 0 gives a constant depth model and HFAC = 1.0 includes 
all of the bathymetry. Any intermediate fraction is also possible if ex- 
perimentation suggests such a case would be of interest. 

To integrate equation (8) the three level Adams - Bashforth method 
is used. That is, writing equations (8) as: 




we use: 

C'(t' + At') = £'(f) + [f g'(t') - \ (10) 

To use this equation we must write g' in finite difference form. This will 
be done term by term starting with g' written as: 

g' = (7 x 0') • V 

+V0’ •?(£,)] + Vx(£.) (ID 
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The first term on the right hand side of equation (11) represents the net 
rate of potential vorticity advected into a unit area. To approximate this 
for the hexagon centered on the point L we write: (all quantities are 



non dimensional and primes have been omitted for simplicity) 



(vx*k) • 7(2^-) - [» L _ N+1 - # L+ 1 >< 



PV +PV 

L-N+l L+l 



, , '( PV L-N + PV L-N+1 ) , ( PV L -1 +PV L-n ) 

^L-N ^L-N+l 2 ^L-l ~ ^L-N 2 

, ,( PV L + N-1 +PV L-1 ), f . . / PV L + N +PV L + N-1 } 

^L+N-l ^L-l' o ^L+N ^L+N-l 2 



+ ( ^L+1 ’ ^L+N } 



(PV + PV ) 
v L+l L+n' 



( 12 ) 



Where: 



and 



PVj = (Cl)(4j) + fj 



J /> 

o 

Cl = a. = —? — 

A DF 



The second term on the right hand side of equation (11) represents the 
diffusion of vorticity horizontally and in finite difference form becomes: 



8V 4 = (C2)T4 + 4 + 4 + 4 

pv S v^^ls l+1 * l _ n+1 s l-N s L -1 



+ 4 +4 

S L+N-1 ^L+N 






(13) 



where 



c2 = 24 = _2K 
3 3A Z F 



the third term on the right hand side of equation (11) represents the 
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dissipation of vorticity through bottom friction and is written as: 





L+N-l 



1 





“ ^L+N-l + ^L-N "^L+N-* 



(14) 




L-N+l 



1 




where: 




The last term on the right hand side of equation (11) is the torque 
added per unit time by wind stress. This is assumed constant in time and 
calculated only once at the beginning of the program using the following 
finite difference form . 



Where G1 is the component of the wind stress in the direction from L to 
L + 1 and G2 is the component of the wind stress in a direction 90 degrees 
to the left of G1 . 

Equations (12) through (15) are used in equation (11) which in turn is 
used in the time integration scheme defined in equation (10). Once the 
vorticity has been obtained for a new time step equation (9) is solved by a 




(15) 




h L-N+l 
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/ 



successive overrelaxation technique. 

The algorithm used is a slight modification of the scheme suggested 
by Winslow (1961). The scheme used is as follows: 

A residual is calculated: 



RES = 7(^ 70) - 4 

^L+l 



= (— ){• 

3 ( h L+1 ~ h L ) 



0 



L-N+l 



0 



L-N 



L-N+l + V (h L-N + V 



(16) 



*L-1 



0 



L+N-l 



L+N 



(h L-l +h L ) (h L+N-l + V (h L+N "L 



+ h~T " 0 l (HFAC l ) 



where: 



HFAC X = (- — ^— ) + ( 



L h L + l +h L 



h L-N+l +h L 



) + (' 



h_ „ T + h_ 
L-N L 



( h +h ) + i h + h } + ( h + h ) 

L-l L L+N-l L L+N L 



The residual is then normalized using the coefficient of 0 in equation 

J_j 

(16) , i.e. , 



= RES 
J_» 



and the relaxation is then given by: 

0 L = 0 L + (R)d0 L (17) 

A number of tests were run to estimate the optimal value for R and over a 
wide range of cases a value of 1.48 was indicated. This value was then 
used throughout the modeling efforts reported here. 

The typical experiment anticipated for the model will be to apply some 
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stress field to a static ocean. The model ocean will then spin-up, or 



develop a circulation pattern that will generally approach a steady state 
providing an appropriate balance of torque is possible and that the 
applied stress field is steady. The time dependent development of the 
steady state circulation and to some extent the final flow pattern will 
depend on the characteristics of the planetary waves that make up the 
transients in the model. For this reason it is of some interest to look 
at the propagation characteristics of these waves within the model and 
to estimate the errors introduced by the finite difference scheme. 

To estimate how the model transients will respond a simple, free, 
linear Rossby wave will be considered. For this case the vorticity 
equation (5) reduces to: 



A wave form of the solution is assumed and near by values used in the 
computational molecule (Fig. 2) taken as follows; 



5 / 1 2 . , f > 

-(-v I/O = -hv • v(^) 



(18) 






i(kx + ly - tot) 





^ (k+</3J) 





qp (-k+i/Ti) 





-ik A 




( 19 ) 
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^L+N-l ” C 

y^k-/^) 

^L+N = * *^L 



Where x and y refer to position in a right handed co-ordinate system 
superimposed on the finite difference grid with the positive x-axis in 
the direction form L to L + 1 . It can also be assumed with no loss of 
generality that the x and y axes are directed east and north respectively- 
This then gives: 



“■ = 0 ; 
ax 



ay" * 



( 20 ) 



To obtain a reasonable estimate of the phase velocity for the Rossby 
waves in the model the time differencing can be assumed exact and the 
depth assumed constant. Using equations (19) and (20) with the finite 
difference forms given in equations (12) and (13), the vorticity equation 
given in (18) becomes, after some simplification. 



• r 1 A _i_ O A* ( k A \ o'! 

-ho [cos kA+ 2 cos (y — ) cos (y- )-3J 



-2M* r . , . , . ,kA . ,,/3£A 

— — — [sin kA+ sin ( — — ) cos (y 



)] 



Which gives 



uj = 



5* A 
2 



sin (kA) + sin (■y ) cos C^y ^ ) 

. cos (kA) + 2 cos (y ) cos ) _3 . 



For small A the trig functions can be expanded using small angle approxi- 
mations . 
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This gives: 



60 = -/?’ 



I(k 3 + k^ 2 )(A 2 ) 

__8 

4 2 2 

,2 2 ,k k f . /a 2 

k + * + Hi H 8 A 



And from this dispersion relation it can easily be seen that the x com- 
ponent of the phase velocity is: 



C 

x 



OJ 

k 



, 1 ,, 2 2 . . 2 , 

1 -- (k +£ )(A ) 



4 2 2 

,2 , 2 ,k k » wa 2. 

+k +( Ti + -8- )(A) 



( 21 ) 



Clearly as ^goes to zero this goes to the exact form 

★ 

R 

C X - - 2 2 

k +ji 

and the error term is second order in A. This then suggests that the 
transient wave solutions that develop in the model will be accurately 
represented to second order. 

The actual fortran program used for the model and an explanation of 
the input data required is given in an appendix. 



Section III - Preliminary Results 

The characteristics of the model described in the last section have 
been investigated in a series of tests designed to check out the numerical 
properties of the solutions. In a number of cases it was also possible to 
illustrate the physical significance that the processes represented might 
have on the circulation in the Arctic Ocean. 

The first test was to check the part of the program that handled the 
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Figure 3. Configuration of the model used to check out the relaxation of the stream 
function. 
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Figure 4. Streamlines of source-sink driven flow with zero vorticity in the interior 
and uniform depth. 
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Figure 5. Streamlines of source-sink driven flow with uniform negative 
vorticity in the interior and uniform depth. 



relaxation of the stream function. A basin shape was chosen that 
roughly covered the deeper portion of the Arctic Ocean (Fig. 3). Repre- 
sentative source-sink distributions were specified as follows: 

03 i 63 i 

2 x 10 m sec flow in across the Chukchi Sea, 7x10 rn sec 

6 3 ^1 

flow in spread out on either side of Franz Josef Land, 1 x 10 m sec 

6 3 “1 

flow out through M'Clure Strait, 1x10 m sec flow out into the 

6 3 _ i 

Lincoln Sea, and 7 x 10 m sec flow out into the East Greenland 
current. With a flat bottom the equations governing the flow in the in- 
terior reduce to: 

V 2 </> = 0 

Solving this for the specified boundary conditions gives the flow indi- 
cated in Figure 4 . 

The next test case introduced some finite vorticity into the above 
equation. For the initial testing this was simply specified as a constant. 
Thus the governing equation becomes: 

V 2 j h = hi 
^ o 

This obviously corresponds to the artifical case where equation 8 has 
gone to a steady-uniform solution for the vorticity. Although this is not 
very realistic the results are interesting and shown in Figure 5 for the 
case where the vorticity in the flow is constant and negative as one might 
expect to get from the stress field developed by the polar easterlies. 
Circulation resembling the Beaufort Gyre develops. Qualitatively the 
results show a striking resemblance to the very viscous solution 
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for ice drift presented by Campbell (1965, Fig. 7c, d). 

The next set of experiments with the model were to check out the 
time dependent vorticity equation (8) in conjunction with relaxation of 
the stream function (Eq. 9). To do this a basin roughly the shape of the 
Arctic was specified. Non-linear terms and bottom friction were not 
included (a = y = 0). A uniform negative stress curl was applied to the 
water that was initially at rest. The Coriolis parameter was assumed 
constant . 

For the first series of runs made on this configuration the depth was 
assumed constant. Under these conditions the model would spin-up 
forming a symmetric clock-wise circulation pattern. The steady state 
solution represented a balance between the torque added by the wind and 
that which was lost through lateral friction. The magnitude of the 
steady state solution and the spin-up time of the model depended on the 
magnitude of the frictional coefficient (8) that was used. It is interesting 
to note that for this series the vorticity equation reduces to: 

^Mv^+J-(7xr) (22) 

dt h 

With a steady wind stress this equation has no wave solutions and the 
appropriate von Neumann type stability condition (Richtmyer and Morton, 
1967) will be of the form: 

St ^ constant » (23) 

(Assuming the Adams-Bashforth scheme is used for time differencing.) 
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o 




Figure 6. Streamlines of flow driven by uniform stress curl in 
an irregular shaped ocean with a central ridge. (The 
ridge runs from top to bottom in the figure.) 
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Where t is the non-dimensional time step in half pendulum days. This 
is relatively weak restriction in that with more or less realistic geo- 
physical parameters time steps of a number of days are possible. This 
was confirmed with the model. ' 

The second series run on this configuration was designed to test the 
models response to variable depth. In this series a ridge of smooth 
profile was placed across the basin. The appropriate form of the 
vorticity equation was 

^-f(vxp) =Pv 2 C + vx(j) (24) 

Under these conditions topographic Rossby waves are possible (Veronis, 
1966) and the von Neumann stability criterion takes on the more restrictive 
form: (Once again assuming an Adams-Bashforth scheme for the time 

differences) 



— < constant 
Af 



1 

2 



(25) 



Where C is the phase speed of the topographic Rossby waves. This 
general behavior was again confirmed by the model. 

In this series the spin-up of the model begins as before, but the 
variable depth acts in many ways analogously to a variable F in a flat 
bottom ocean and western boundary type currents can develop. These 
stronger boundary currents develop not on the western side of the ocean 
basin as in Munk's model (Munk, 1950) but rather to the left of an observer 
looking from deep to shallow water. Figure 6 shows a typical solution 
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for this series. The ridge is seen to divide the flow into two cells 
with regions of intensified currents in each half. 

During the initial spin-up for this series transient wave patterns 
were seen to propagate along the ridge. It is quite likely that the real 
time dependent stress fields applied to the Arctic ocean will result in 
transients of this form being propagated along the Lomonosov Ridge 
(Fig. 1). 

An additional series of tests were run to investigate the interaction 
of source-sink terms with the bathymetry and the formulation used to 
represent bottom friction. For this series the stress applied at the 
surface was zero. The source-sink distribution was simply to have flow 
in through the Bering Straits and out through the East Greenland area. 

The bathymetry was again the smooth ridge described in the previous 
series of tests . 

The first test in this series assumed oc =3 = y = 0. Thus there was 
no friction, or non-linear terms. The appropriate form of the vorticity 
equation that was: 

- (vx i|)k) * V (j^) = 0 (26) 

This equation shows the interesting result that the only steady flow 
possible is when the stream lines are parallel to lines of constant f/h. 
Thus the steady-state flow can not cross the ridge and any source-sink 
distribution that demands cross ridge flow will not lead to a steady 
circulation. The numerical model shows these characteristics with 
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circulation continuing to increase with time. A strong negative (clock- 
wise circulation builds up on the in flow side of the ridge and a 
corresponding positive circulation develops on the outflow side of the 
ridge. This is similar to the effect one might get if the water simply 
piled up on the in-flow side and drained out on the out-flow side 
satisfying the overall continuity requirements for the ocean but with 
minimum cross-ridge flow. 

This has an important implication for the circulation in the Arctic. 

If the source-sink component of the flow is to be baratropic and geostrophic 
then flow can not cross the Lomonosov ridge. This means that continuity 
must be satisfied separately for both basins in the Arctic. There will 
be a strong tendency for flow that enters on the Canadian side of the Arctic 
to exit on the same side. Field data from the Arctic suggest that this may 
in fact be the case (Coachman, 1970). Thus one might expect that flow 
through the Bering Straits and the Canadian Archipelago should be strongly 
correllated . 

The next runs in this series included the non-linear terms and 
lateral friction which gave the following governing vorticity equation: 



For this case both frictional and non-linear boundary layers are present 
with cross ridge flow taking place and a steady-state solution is possible. 
For all reasonable values of a and 3 the interior of the flow is still essen- 
tially geostrophic and the cross ridge flow takes place in boundary layers 




(27) 



30 



Figure 7. Streamlines of flow driven by a simple source-sink 

distribution for an irregular shaped ocean with a central 
ridge. (The ridge runs from the top to the bottom in the 
figure . 
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where the ridge intersects the sides of the basin. (Fig. 7) 

A final case run in this series included the bottom friction terms. It 
can be anticipated that bottom friction would be effective in causing 
cross bathymetry flow. In particular the 

r’* • 7< h ) 

term will have a tendency to turn flow towards shallow water. This 
is in fact what the model showed. Steady state flows resembled figure 7 
for very small values of y and showed a smooth transition to stream line 
patterns resembling simple hydraulic flow (similar to Fig. 4) as the value 
of y was increased. 

In the deeper interior of the Arctic Ocean the flow must be essentially 
geostrophic, but it seems likely that frictional effects must be significant 
over large portions of the Siberian shelf. In the vicinity of the Lomonosov 
Ridge both friction and non-linear effects may be significant. 

Section IV - Future Model Plans 

The next stages in the development and use of this model are well 
under way at the time of this writing. 

The first extension is to use the model with greatly increased reso- 
lution and the actual bathymetric variation of the Arctic basins. With 
this configuration realistic source-sink distributions are investigated. 

A wind stress field obtained from Felezenbaum's pressure data 
(Felezenbaum , 1958) is used to drive the flow and a number of runs are 
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anticipated to investigate the effects of variations in a, 3 and y . 

The next series of tests will use the same basic model configuration 

but the stress field will be obtained by applying Felezenbaum's pressure 

* 

data to Campbell and Rasmussen's (1971) ice model and then using the 
stress field from the bottom of the ice to drive the ocean model. Again a 
series of runs are anticipated to investigate variations in parameters. 

Another series of runs is planned for the model that will have even 
higher spacial resolution and definitions of the bathymetry. In this 
series the ocean will be driven by a time dependent wind stress field that 
is made up of random frequency components. The dependent variables at 
a number of locations will be spectral analyzed with the intension of 
obtaining normal mode frequencies for the Arctic basins. 

The next stage in the development of model exploration for the Arctic 
will be the formulation of a baroclinic model. The first set of experiments 
anticipated for this formulation will be the investigation of the dynamics 

associated with the circulation of the Atlantic Layer (Coachman and Barnes, 

/ 

1963). 



\ 
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Appendix - Use of Computer Program 



The ocean model program has been written in FORTRAN - IV and a 

listing of the program is included in the end of this appendix. To use 

the program an appropriate data deck must be placed after the program. 

The composition of the data deck must be as follows: 

1) N, M, NB FORMAT (313) - this is a single card where N and M are 
the dimensions of the grid system (ref. Fig. 2) and NB is the number 
of boundary points. 

2) IA FORMAT (2013) - this is a grid number list of the first interior 
point of each line down the left hand side of the model. 

3) JA FORMAT (2013) - this is a grid number list of the last interior 
point of each line down the right hand side of the model. 

4) LV, BLANK FORMAT (21A1) - this is a single card with the code to be 
used in the graphical out put of the program. (Ref. subroutine GRAOUT) 

5) HM, FAC FORMAT (2F6.1) - this is a single card where HM is the 
reference depth and FAC is the fraction of the bathymetry that is 
desired in the model. 

6) H FORMAT (12F6.4) - this is a set of cards that has the normalized 
depth variation for each grid point. 

7) F FORMAT (12F6.4) - this is a set of cards that has the normalized 
values of the Coriolis parameter for each grid point. 
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8) I FORMAT (13) - this is a single card with the exponential part of 
the values to be used for the wind stress. 

9) GI FORMAT (12F5.3) - this is a set of cards with the significant 
figures part of the wind stress in the direction of the positive axis, 
i.e. from point L to L + 1 (Ref. Fig. 2). 

10) G2 FORMAT (12F5.3) - this is a set of cards with the significant 
figures part of the wind stress in the direction 90° to the left of the 
positive axis . 

11) IBC FORMAT (6(611, 13)) - this is a set of cards that contains the 
information required to calculate boundary values of vorticity. For 
each broundary point a six digit code and the grid number are given. 
The six digit code refers to the neighboring points in a counter clock- 
wise order starting with L + 1. The following code is used: 

0 - exterior point, 1 - interior point (free slip B.C.), 2 - interior 
point (no slip B.C), 3 - boundary point along streamline, 4 - boundary 
point across a streamline. Boundary points at sources are not in- 
cluded in the IBC list. 

12) Cl, C2, C3 , DT, TOUT, TMAX, R, CON, 1PUNCH FOPMAT(6E10 .3 , 
F4.2, E10.3, 12) - this is a single card that contains the run para- 
meters. Cl, C2 and C3 are defined in section two of this report. DT 
is the time step. TOUT is the time interval at which output of the 
dependent variables is desired. TMAX is the maximum time the model 
is to run. R is the SOR parameter to be used for the relaxation of the 
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stream function. CON is the absolute convergence limit to be used 



for the relaxation of the stream function. IPUNCH is a flag that 
should be non-blank if the final values for the dependent variables 
are required as a source deck for future runs. 

13) S FORMAT (6E12.5) - this is a set of cards that contain the initial 
values for the stream function at each grid point. 

14) V FORMAT(6E12 . 5) - this is a set of cards that contain the initial 
values for the vorticity at each grid point. 

15) G1 FORMAT(6E12 . 6) - this is a set of cards that contain the initial 
values of the local time rate of chance of the vorticity at time -(-DT). 
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